
rm(list=ls())

library(tidyverse)
library(tidycensus)
library(stringr)
library(stargazer)
library(lubridate)
library(viridis)
library(sf)
library(sp)
library(geosphere)
library(writexl)
census_api_key("14cdfe8a1ca6db7607349a73135e1329c2397813")

Boundaries_Path <- "../Data/Boundaries/"
# Mapping
library(raster) # Shape File Manipulation
library(rgeos) # Geospatial
library(rgdal) # Geospatial
library(ggmap) # Google Maps Layers

archive <- Sys.Date() %>%
  ymd() %>%
  str_replace_all(., "-", "")

theme_opts <- list(theme(panel.grid.minor = element_blank(),
                         panel.grid.major = element_blank(),
                         panel.background = element_blank(),
                         plot.background = element_rect(fill=NA),
                         panel.border = element_blank(),
                         axis.line=element_blank(),
                         axis.text.x=element_blank(),
                         axis.text.y=element_blank(),
                         axis.ticks=element_blank(),
                         axis.title.x=element_blank(),
                         axis.title.y=element_blank(),
                         plot.title = element_text(size=22)))


## Projection definition
ACS_projection <- st_crs("+proj=longlat +datum=NAD83 +no_defs")

# Census ACS data
zcta_shapefile <- get_acs(geography = "ZCTA",
                              variables = "B19013_001", 
                              year = 2013, 
                              output = "wide",
                              state = "CA", 
                              geometry = TRUE) %>%
  st_transform(., ACS_projection)

# Collect Census data
bg_2013 <- get_acs(geography = "block group", 
                   variables = c(income = "B19013_001", 
                                 population = "B01001_001"),
                   output = "wide",
                   state = "CA", 
                   geometry = TRUE) %>%
  st_transform(3488)


service_areas_sf <- st_read("../Data/Utility_ShapeFiles/Electric_Utility_Service_Areas") %>%
  st_transform(., ACS_projection)
#service_areas_df.no_sf <- as.data.frame(service_areas_sf) # in case we want the data in dataframe, not sf. 
service_areas_sf_exp <- service_areas_sf
st_geometry(service_areas_sf_exp) <- NULL
write_csv(service_areas_sf_exp, "../Data/Boundaries/service_areas_sf.csv")
rm(service_areas_sf_exp)
geom <- st_geometry(service_areas_sf)

# Add area to CBG data
bg_2013 <- bg_2013 %>% 
  mutate(area = st_area(.),
         area = as.numeric(area))

bg_2013_exp <- bg_2013
st_geometry(bg_2013_exp) <- NULL
write_csv(bg_2013_exp, "../Data/Boundaries/bg_2013_exp.csv")

med_income_bg_2013_exp <- bg_2013
st_geometry(med_income_bg_2013_exp) <- NULL
write_csv(med_income_bg_2013_exp, "../Data/Boundaries/med_income_bg_2013.csv")
rm(med_income_bg_2013_exp)

############################################################################################
## create datasets on overlaps and distances, 
############################################################################################

##but first st_transform
med_income_bg_2013 <- bg_2013
geom <- st_transform(geom,3488) 
med_income_bg_2013 <- st_transform(med_income_bg_2013,3488)
zcta_shapefile <- st_transform(zcta_shapefile, 3488)
med_income_bg_2013_geo <- st_transform(med_income_bg_2013, 3488)

# Map zctas to cbgs
join_bg_alt_stations <- st_join(med_income_bg_2013, zcta_shapefile)%>% data.frame() %>% write.csv("../Data/Boundaries/join_ZCTA_bg.csv")

# create centroids of the 2015 Census tracts, distances between tracts and utility boundaries
bg_2013_centroid <- st_geometry(med_income_bg_2013) %>% st_centroid()
bg_2013_area <- st_geometry(med_income_bg_2013) %>% st_area() %>% data.frame() %>% write_csv(bg_2013_area,"../Data/Boundaries/bg_2013_area.csv")

# distances between tracts and utility boundaries
distance2 <- bg_2013_centroid %>% 
   st_distance(geom)
data2 <- data.frame(distance2) %>%  write_csv("../Data/Boundaries/distance_data_2013_bg.csv")


####
## THIS TAKES A LONG TIME. ONLY RUN IF NECESSARY
####

seq <-  seq(1,12)
for (val in seq) {
  lower = (val-1) * 2000 + 1
  if (val != 12){
    upper = val * 2000
  } 
  else {
    upper = 23212
  }
  bg_2013_centroid_sub <- bg_2013_centroid[lower:upper]
  name <- paste("../ResultsOut/Boundaries/distance_bg_2013_sub",val,".csv", sep ="")
  distance_bg <- st_distance(bg_2013_centroid, bg_2013_centroid_sub) %>% data.frame() %>% write_csv(name)
  rm(distance_bg, bg_2013_centroid_sub)
}

####
# adjacency for utility boundaries
adjacent_utilities <- st_intersects(service_areas_sf, service_areas_sf)%>% data.frame()%>% write_csv("../ResultsOut/Boundaries/adjacent_utilities.csv")


overlap_boundaries_proper <- st_contains_properly(geom, med_income_bg_2013) %>% data.frame %>% 
  write_csv("../ResultsOut/Boundaries/overlap_boundaries_proper_bg.csv")
overlap_boundaries_proper_new <- st_join(bg_2013, service_areas_sf)%>% data.frame() %>% write_xlsx("../ResultsOut/Boundaries/overlap_boundaries_proper_new.xlsx")


